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A one-dimensional model of inertial pumping is introduced and solved. The pump is driven by a 
high-pressure vapor bubble generated by a microheater positioned asymmetrically in a microchannel. 
The bubble is approximated as a short-term impulse delivered to the two fluidic columns inside the 
channel. Fluid dynamics is described by a Newton-like equation with a variable mass but without 
the mass derivative term. Because of smaller inertia the short column refills the channel faster and 
accumulates a larger mechanical momentum. After bubble collapse the total fluid momentum is non- 
zero resulting in a net flow. Two different versions of the model are analyzed in detail, analytically 
and numerically. In the "symmetrical" model, the pressure at the channel-reservoir connection 
plane is assumed constant whereas in the "asymmetrical" model it is reduced by a Bernoulli term. 
For low and intermediate vapor bubble pressures, both models predict the existence of an optimal 
microheater location. The predicted net flow in the asymmetrical model is smaller by a factor of 
about two. For unphysically large vapor pressures, the asymmetrical model predicts saturation of 
the effect, while in the symmetrical model net flow increases indefinitely. Pumping is reduced by a 
non-zero viscosity but to a different degree depending on the microheater location. 

PACS numbers: 47.60. Dx, 47.61. Jd 



I. INTRODUCTION 

In a recent paper [l[ , we reported numerical modeling 
and experimental demonstration of inertial pumping. In 
this effect fluid is confined within a thin channel 

connecting two large reservoirs and pumped by a high 
pressure vapor bubble that repeatedly expands and col- 
lapses. The bubbles can be created by a number of actua- 
tion methods including thermal resistors [l|, |3| , electrical 
currents passing through the fluid Jl, Q , localized laser 
pulses p-|3| , and acoustic actuation Q . If the bubbles are 
generated away from the geometric center of the channel, 
the dynamics of the two fluidic columns are asymmetric, 
which results in a net flow from the short toward the long 
arm of the channel. 

The pumping action has been attributed to unequal 
inertia of the two fluidic columns inside the channel [l, d, 
0, • While several numerical studies have been done in 
realistic three-dimensional geometries [11, d, , the main 
effect can be captured within a simple one-dimensional 
model derived from momentum balance [H-ll, @1- The 
model has been used to prove nonzero net flow Q and 
to illustrate the main qualitative features of the pump 
through numerical analysis 

The model is nonlinear and in general not solvable ana- 
lytically. However, for realistic values of pressure (several 
atmospheres), channel width (tens of microns), and fluid 
properties (density, viscosity, and surface tension of wa- 
ter), the dynamics is dominated by the balance between 
inertia and pressure forces [l[ . If only those two factors 
are left in the dynamic equation, the latter simplifies so 
much that an analytical solution becomes possible. The 
goal of this paper is to present this solution. It provides 



* |pavel.kornilovich@hp.com| 



new insight into the mechanism of inertial pumping and 
becomes a useful tool for studying this effect. 



II. THE MODEL 

Consider the geometry of Fig. [TJa). A thin channel 
of cross-sectional area A is connected to two much wider 
reservoirs. The reservoirs and the channel are filled with 
incompressible fiuid of density p. The pressures in the 
reservoirs far from the channel are pit and p2b, respec- 
tively. In the most common case, both bulk pressures 
are equal to the atmospheric pressure po but more gen- 
eral cases may be considered. Inside the channel is a 
resistive microheater that locally boils the fluid and cre- 
ates high-pressure vapor bubbles. Each bubble expands 
and collapses resulting in a non-zero net flow through the 
channel, see Fig. [Ijb)-(g). The microheater must be po- 
sitioned asymmetrically relative to the channel's ends to 
achieve pumping. 

The simplest approach to flow dynamics is to com- 
pletely neglect transverse motion of the fluid in the chan- 
nel, three-dimensional features of the bubble shape, and 
curvature of the vapor-fluid interface. These assump- 
tions are approximate; however, they do not change 
the qualitative picture of inertial pumping. Through 
full three-dimensional Computational Fluid Dynamics 
(CFD) modeling we showed that pumping is most effi- 
cient when the bubble occupies the entire cross-section of 
the channel pushing the fiuid primarily along the chan- 
nel's axis Under these assumptions, the channel is 
approximated by a one-dimensional line and a vapor-fiuid 
interface by a single point on this line. Thus a single bub- 
ble is fully described by the time-dependent positions of 
two interfaces xi and X2. The case of more than one 
bubble can be formulated in the same manner. The goal 
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FIG. 1. (a) Schematic geometry of an inertial micropump. 
(b)-(g) Phases of the expansion-collapse cycle, (b) The start- 
ing state: the fluid is at rest, (c) Microheater (the black 
square) creates a high-pressure vapor bubble. A positive pres- 
sure difference pushes the fluid out of the channel, (d) The 
vapor bubble pressure quickly drops below atmospheric, and 
the fluid decelerates under a negative pressure difference while 
moving by inertia, (e) The short arm turns around while the 
long arm is still moving out of the channel, (f ) The two fluidic 
columns collide at a point shifted from the expansion starting 
point: the primary pumping effect, (g) Since the shorter arm 
has a larger momentum at collision, the total post-collapse 
momentum is non-zero: the secondary pumping effect. 



of the present paper is to analyze the single bubble case. 



A. Momentum balance and forces 

Consider the left column of fluid (see Fig. [T]). Its 
length is Xi{t), the velocity is xi, and total momentum 
is Qi{t) = pAxi xi. In a time interval dt the momentum 
changes because of two factors: (i) force Fi{t) acting on 
the fluid, and (ii) momentum lost to or supplied by the 
left reservoir. During the expansion phase xi < 0, and 
the column loses a negative momentum pAx\ dt to the 
reservoir. This element must be added to momentum 
balance as a positive increment: 

Qi{t + dt) = Qi{t) + pAx\dt + Fi{t)dt . (1) 

During the collapse phase xi > 0, and the column ac- 
quires a positive momentum pAxi dt from the reservoir 
and the same equation applies. Substituting Qi one de- 



rives the dynamic equation 

pAxiii^Fiit). (2) 

This is a variable-mass Newton equation but without 
the mass-derivative term: the latter has been dissipated 
by the reservoir. A similar momentum balance analysis 
yields for the right fluidic column 

pA{L-x2)!i2=F2{t) , (3) 

where L is the total channel length. In this case, the 
right reservoir absorbs or supplies the extra mechanical 
momentum. 

The major forces at play during the expansion-collapse 
cycle are: (i) pressure-difference forces, (ii) viscous forces, 
and (iii) surface tension forces. The surface tension forces 
are about two orders of magnitude smaller than the other 
two for typical microfluidic conditions [ij and can be 
safely neglected. The viscous force is typically smaller 
than the pressure force but is not always negligible. The 
simplest way to account for the viscous force is to assume 
that it is proportional to velocity and also to the length of 
the column because it acts along the fluid-surface bound- 
ary. The resulting one-dimensional model reads 

Apxixi + KXiXi ~ {pi - Pv)A , (4) 

Ap [L — X2) X2 + k{L — X2) X2 = {pv — P2)A . (5) 

Here k characterizes the relative strength of the viscous 
force. Under further assumptions, k can be linked to 
other parameters of the system. For example, within the 
Poiseuille flow 

K = Sirr/ , (6) 

where rj is the bulk viscosity of the fluid. This relation- 
ship is derived in the Appendix. In this work k will be 
treated as an independent phenomenological parameter. 

B. Impulse bubbles and boundary conditions 

The next step is to specify the pressures in Eqs. (|4])-([5|). 
There, Pv is the vapor pressure inside the vapor bubble 
while pi^2 are the fluid pressures inside the reservoirs near 
the channel's ends, see Fig.[IJa). All pressures are time- 
dependent and require separate considerations. Let us 
begin with the vapor bubble. 

A typical pumping event starts with a microheater 
boiling the surrounding fluid within a fraction of a mi- 
crosecond. This creates a thin layer of superheated vapor 
with pressure of about 8-10 atm This high-pressure 
phase lasts about one microsecond. Then the pressure 
quickly drops to sub-atmospheric levels because of bub- 
ble expansion and heat transfer losses. The residual bub- 
ble pressure (the saturation pressure of the vapor at the 
ambient temperature) is about = 0-3 atm. Since the 
entire pumping cycle lasts dozens of microseconds, the 
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initial high-pressure phase can be approximated as in- 
stantaneous impact on the fluid, as a resuh of which the 
fluid acquires mechanical momentum go This leads 
to the following impulse bubble model: 

PvW = f^W+Pvr. (7) 

In this approximation, the bubble strength is character- 
ized by the initial fluid momentum qg. Since the impulse 
has negligible time duration, qo drops out from the model 
equations and enters the dynamics only via the initial 
conditions. This greatly simplifles analysis. If a detailed 
physical model of the vapor bubble is known, qo can be 
estimated from Eq. ([7]) by integrating the vapor pressure 
over the bubble's lifetime. In this work, qo will be treated 
as an independent phenomenological parameter. 

Consider now the pressures at the channel's ends pi,2- 
In general, they are different from the respective bulk 
pressures pu and p2b- The fluid dynamics in the reser- 
voirs is complex. First, the formation of a jet during 
outward flow and the formation of a sink during inward 
flows results in asymmetry. Second, the pumping event is 
inherently transient. The fluid expelled in the reservoirs 
forms three-dimensional vortices that persist long after 
the bubble collapse [3| . One consequence of the vortices 
is non-homogeneous inward flow: the fluid is filling up 
the channel along the walls while emptying the channel 
near its axis at the same time. All of this makes the 
task of finding a representative one-dimensional bound- 
ary condition non-trivial. However, a detailed analysis of 
the reservoir flow is beyond the scope of the present work, 
and is deferred to a separate study. The main goal here 
is to elucidate the physics behind inertial pumping by us- 
ing simplified one-dimensional dynamics as the analysis 
tool. From this perspective, replacing the complex reser- 
voir dynamics with a deterministic boundary condition 
at the end of the channel is acceptable, as long as it is 
physically justified and consistent with net pumping. 

The easiest approach is to simply neglect all the above- 
mentioned effects and set the end pressures equal to the 
bulk values: pi = pib and p2 = P2b- This choice will be 
called the symmetrical model because it docs not distin- 
guish between outflow and inflow. The main advantage 
of the symmetrical model is simplicity; it is the easiest 
way to derive the main qualitative features of inertial 
pumping. 

A different model was proposed by Yuan and Pros- 
peretti who tried to account for the asymmetry be- 
tween expansion and collapse. During bubble expansion, 
the fluid is injected into the reservoir as a jet that sep- 
arates from the rest of the fluid. The pressure near the 
entry point is close to the bulk pressure of the reservoir, 
Pi = Pib- During collapse, however, part of the pres- 
sure difference is expended on accelerating the fluid in- 
side the reservoir. This is a "sink" flow According 
to the Bernoulli equation, pi = pib — ^pif- As mentioned 
above, vortex formation near channel ends disturbs the 
sink flow, and it is unclear to what extent the Bernoulli 



correction can be applied to the situation at hand. In the 
absence of detailed analysis of the pressure distribution, 
the proposal of Yuan and Prosperetti will be adopted 
here as the second reasonable choice of the boundary 
condition. Hereafter, it will be referred to as the asym- 
metrical model. Compared with the symmetrical model, 
it provides additional insight into the inertial pumping 
mechanism, but is more complex mathematically. 

It is of benefit to keep both options available for 
analysis. To this end, wc introduce the discrete index 
m = {0, 1} that will distinguish between the models: 
m = corresponds to the symmetrical model and m = 1 
to the asymmetrical model. With such notation the pres- 
sure boundary conditions can be written as 

Pi =Pib- — pilH{xi) , (8) 

P2^P2b-^ P±lH{~X2) ■ (9) 

Here H{z) is the Heaviside step function: H{z > 0) = 1, 
and H[z < 0) = 0. 



C. Dimensionless pump equations 

The number of model parameters in Eqs. (Il))-® can be 
reduced by transforming the equations to a dimension- 
less form. In a similar procedure in Ref. [l|, the channel 
diameter was taken as the unit of length while the dura- 
tion of the high-pressure phase as the unit of time. In the 
present case, both these parameters are zero, so a differ- 
ent pair of units is required. The choice of unit length 
is obvious: the total channel length L is the only model 
parameter of that dimensionality. To choose a unit time 
one can proceed as follows. In the most typical case, bulk 
reservoir pressure is equal to the atmospheric pressure po ■ 
The two fluid columns will be evolving under negative 
pressure difference — (po ~ Pvr)- That sets a character- 
istic fluid velocity as (po — Pvr)/p and a characteristic 
time as L\J pj (po — Pvr)- The characteristic time thus de- 
flned is of the order of the total bubble lifetime and the 
characteristic velocity is approximately mean fluid veloc- 
ity over the same period. Accordingly, one introduces 
dimensionless time r and interface positions S,i 

= ^ . (11) 

Upon substitution of Eqs. ((7|)- (|TT|) . the pump equations 
Q-® become 

CiCi +yefi?(Ci)+/56Ci =71, (12) 
(1 - 6) - Y^'iHi^Q + (1 - 6) ^2 = -72,(13) 
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where < £,i{t) < ^2(,t) < 1, the primes denote differen- 
tiation with respect to r, and 



P 



71,2 



pAy po- pvr ' 
Plb,2fc - Pvr 
Pa - Pvr 



(14) 

(15) 



The equations of motion contain one discrete model pa- 
rameter m = {0, 1} and three positive continuous dimen- 
sionless parameters: friction /3 and pressures 71, 72. In 
the most common case of bulk pressures equal to the at- 
mospheric pressure, 71 = 72 = 1, only one continuous 
parameter /3 remains. 

Consider now initial conditions. The expansion phase 
starts with a zero-size bubble located at xq, hence 



ei(o) = 6(o) = ^ = eo. 



(16) 



The microheatcr location ^0 is a critical parameter of 
the inertial pumping effect. For a symmetrically-placed 
microheater, ^0 = 1/2, net flow must be zero. Net flow 
increases as ^0 deviates from 1/2 but in a non-monotonic 
manner. This behavior will be analyzed in the following 
Sections. 

Another important parameter is the bubble strength. 
Within the impulse-bubble approximation, at t = the 
two columns of fluid acquire momentum qq. Converting 
this to initial velocities, one obtains wi(0) = —qo/pAxo 
and ^2(0) — qo/pA{L — xq). In dimensionless units 



a 



1-^0 



where 



pAL y Po- Pvr 



(17) 



(18) 



The initial conditions, Eq. (|17p . assume that expansion 
begins with a state of rest. One may consider a more 
general case of fluid already moving as a whole at r = 
0. (It occurs for example with high-frequency repeated 
pumping or when pumping against an imposed pressure 
gradient.) This case is not included in the present work 
and is deferred to a separate study. 

To summarize, bubble kinematics is completely defined 
by the positions of vapor-fluid interfaces satisfying < 
?i < < 1- Bubble dynamics is governed by equations 
of motion ((T^ - ([T^ and initial conditions ([T5|) and ([TT]) . 
There are five continuous dimensionless parameters: a, 
/3j 71, 72, £,0i and one discrete model index to. Typical 
parameter values are given in Table ID 



D. Post-collapse phase 

The collapse phase ends with the two fluidic columns 
colliding at time tc at a point Xc which is different from 



Parameter 


Meaning 


Typical values 


a 


Bubble strength 


0.1 - 1.0 


P 


Friction 


1 - 10 


71.2 


Pressure difference 


1 - 2 


Co 


Asymmetry 


(0,1) 



TABLE I. Physically relevant ranges of the dimensionless 
model parameters. The estimates have been obtained for p = 
10^ kg/m^, L = 200 pra, A = 20^ pm^ , Po = 1 atm and 
Pvr = 0.3 atm P]. Using high-pressure phase duration of 1 ps 
and pressure of 8 atm the initial impulse is estimated as 
qa/A ^ 0.7 kg/(m-s). Then Eq. ^ yields a ^ 0.42, which 
is a typical value of the bubble strength. For a fluid with 
bulk viscosity rj = 1.3 mPa-s, Eqs. (fTi)) and ^ yield /3 ^ 2.0, 
which is adopted as a typical value of the friction coefficient. 



the starting point of expansion xq- Shift Ax = Xc — Xq 
constitutes the primary pumping effect. In addition, mo- 
menta of the two flows are in general different. (This is 
despite the fact that motion started with equal momenta 
qo- Reasons for the eventual inequality are analyzed in 
the next Section.) By momentum conservation, this im- 
plies that after collapse the fluid will continue to flow as 
a whole (i.e. with total length L and mass pAL) until it 
is brought to a complete stop by pressure difference and 
viscous forces. Velocity Vc in the beginning of this phase 
will be referred to as post-collapse velocity and net flow 
during the post-collapse phase as the secondary pumping 
effect. Flow kinematics can still be described by time 
evolution of collapse point x{t). Acting forces no longer 
include the bubble vapor pressure but only end-channel 
pressures pi^2- The equation of motion reads 



pAL X + kL X = (jpi — P2)A 



(19) 



with initial conditions x{tc) = Xc and x(tc) = Vc- Sub- 
stituting the pressure model, Eqs. (H))-® and converting 
to the dimensionless units, one obtains 



r + ^C"sgn(e') + /3e'=7i- 72, (20) 



with ^(tc) = and ^'(tc) = r]c. Parameters in Eq. ([20)1 
have the same meaning as during the expansion-collapse 
cycle. 



III. 



INERTIAL PUMPING EFFECT 



Quite generally, inertial pumping happens because the 
shorter fluidic column with its smaller mass turns around 
faster than the longer one under the same pressure differ- 
ence. The shorter column returns to the starting point 
earlier, and has an extra time before collision to keep 
moving in the long arm's direction. This results in the 
primary effect. For the same reason, the shorter arm has 
more time to accelerate during collapse. Consequently, 
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FIG. 2. Dawson integral, Eq ((25}, its derivative, and the 
function -F{z). Asymptotic expansions are F{z <^ 1) ~ 



z- ^z^ 

^ 3^ ' 15 



+ ±z'' + ... and F(2 > 1) ^ ^ 



+ 



it arrives at collision with a larger momentum than the 
longer arm, loading to the secondary effect. Qualitative 
features of inertial pumping can be understood by an- 
alyzing return times and interface trajectories near the 
collision point. In this Section we study the basic inviscid 
pumping event with /3 = and 7i = 72 = 1. 



A. Expansion phase 

It is sufficient to follow the dynamics of only one vapor- 
fluid interface. Math is slightly more transparent for ^1. 
The initial value problem reads [cf. Eq. ([T^ for < 0] 

iiC^ = 1 

ei(0)=eo . (21) 

Ci(o) = -^ 

Integrating once, one obtains a first integral: 

^ef-ine.^c = i(|)^ineo. (22) 

The interface undergoes potential motion in a profile 
— In^, with C being analogous to a "total energy." At a 
turning point = 0, from where the coordinate of the 
turning point is 



Cit = Co e 



(23) 



Integrating second time from the start to the turning 
point, the total expansion time is given by 



Tie 



= 



e 



dz . 



(24) 



Introducing Dawson integral (Fig. [2]) 

Jo 

the expansion time can be written as 

1 



Tie 



F 



(25) 



(26) 



The expansion time of the right interface is obtained from 
here by substitution ^ 1 ^ Co^ 



1 



T2e 



V2{l-io) 



F 



V2{1 - Co) 



(27) 



A crucial observation now is that the function -F{z) 
monotonically decreases with its argument, as shown in 
Fig. [2l Clearly, expansion times are the same for the 
symmetrical microheater position = 1/2- However, for 
asymmetric positions (and the same bubble strength) the 
expansion times will be different: Tie < T2e for < 1/2, 
and Tie > 'T2e for ^0 > 1/2- Thus, due to its smaller mass, 
the shorter arm slows down faster and reaches the turn- 
around point earlier than the longer one. This forms the 
basis for inertial effects of both types. 



B. Collapse phase. Symmetrical model (m — 0) 

In the symmetrical model collapse dynamics is an in- 
verse of expansion dynamics. Acceleration can be de- 
scribed as sliding along the same potential profile — In ^ 
with a zero starting velocity. The return times (i.e. times 
after which both interfaces return to the starting point 
Co) are simply twice the expansion times, Tir,2r = 2Tie,2e- 
Both the primary and secondary pumping effects can be 
derived analytically for weak and strong asymmetries. 
These cases are considered first. After that a complete 
numerical solution of the pump equations is presented. 



1. Weak asymmetry 

Under weak asymmetry, the microheater is close to the 
channel's midpoint: 



Co = ^ + V' 



ip < 1 



(28) 



We are interested in the pumping effects that are of first 
order in i/j. Substituting in the return times and ex- 
panding to first order, one obtains 



Tir,2r = a±bip + 0{ljj'^) 
a = V2F{V2a) , 



b ^ 2^/^F{V2a) ~ 4aF'{V2a) = 
= 2^/2(1 + 4:a^)F{V2a) - 4a 



(29) 
(30) 

(31) 
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The last line follows from the identity F'{z) = l — 2zF{z). 
The limit values for weak and strong bubbles are 

a(a<l)«2a- — , a(a > 1) w — , (32) 
3 2a 



b{a < 1) 



6(a> 1) 



(33) 



from where asymptotes of the return times can be de- 
duced. Of primary interest are the time and coordinate 
of the collision point. They can be found by linearizing 
interface trajectories near the starting point just prior to 
collision. At r = Tir, the left interface is at = and 
has velocity S,[ = a/(i + V')- Accordingly, a linearized 
trajectory is described by the equation 



1 



a 



(t — a — 6'0) . 



(34) 



Similarly, the right interface's trajectory linearized aro- 
und the collision time is 

6(r) = Q + 4^ ^ (r-a + b^). (35) 

Equating ^i{tc) = ^2iTc) yields the collision time and 
position in the first order in ip: 

T^ = a + Q-^ + 0{tP^) , (36) 
^c= (l+A ~2abij. (37) 



The primary pumping effect is the displacement of the 
collapse point relative to the starting point of expansion. 
According to the last equation, the primary effect is given 
by 



a > 1 



(38) 



Note that the primary effect is a sharp function of the 
bubble strength (cx a*) for weak bubbles. 

The secondary pumping effect originates from the im- 
balance of mechanical momenta at collision. The veloci- 
ties can also be found from linearizing them around the 
collision point. Accelerations follow from the original dy- 
namic equations. The linearized velocities are 



— 7 + T— -(^c-ri^)^^— — , (39) 

—a 1 , , —a — bTp 

-{T,~T2r)^—, —. (40) 



4> 



The total momentum after collapse to the first order 
in if) is 

9c = • + (1 - Cc) • = -2^1 + 4a')V' 



3 

-16 aip . 



a > 1 



(41) 



Since the total "mass" after collapse is equal to 1, the 
post-collapse velocity is given by the same expression. 



The first term in Eq. (|41|) originates from additional ve- 
locity that the short arm acquires between the return 
time and collision time. The second term comes from 
an additional mass of the short arm acquires during the 
same time interval. Depending on bubble strength a, 
either the first or the second factor dominates the sec- 
ondary effect. 



2. Strong asymmetry 

Consider now the opposite case of strong asymmetry 
when the microheater is very close to one edge of the 
channel, say, to the left one. Then the starting point 
Co ^ 1 is a small parameter. Of interest are the primary 
and secondary pumping effects in the lowest order of Co- 

Let us start with the return times. Clearly, it takes the 
short arm a much smaller time to return than the long 
arm. Referring to Eq. (j26p . and using the large- argument 
expansion F{z ^ 1) = (2z)^^ + 0(z~'^), one obtains 

ri,.(eo « 1) = 2rie(Co « 1) = ^ . (42) 

a 

Thus the short arm's return time is second order in Co- 
One power comes from a short travel distance and an- 
other power from a large starting velocity. In contrast, 
the return time of the long arm is 0(1). 

The collision time and position can again be found 
from linearized trajectories around the collision point. At 
r = Tir, the short arm has velocity a/Co acceleration 
1/Co- The trajectory is 

ei(r) = Co + f (r - T,.r) + -^(r - T,,f + ... (43) 
Co ^Co 

The right arm moves comparatively slow; therefore it is 
sufficient to expand its trajectory around the starting 
point and time. It reads: 



6(r) =Co + 



1 



1 - Co 2(1 - Co) 



(44) 



By setting £,i{tc) = C2(7'c) after some algebra one obtains 
the time and position of the collision point 



Tc = Tir{l +Co) 
Cc = Co + 2Co' + 



= ^(1 

a 



-Co) + O(^o'), (45) 
(46) 

(47) 



Thus the primary pumping effect is 

AC = Cc - eo = 2Co' , 

which is quadratic in Co- It is shown in Fig. [3] by the 
dashed line. 

Velocities needed for the secondary effect can be found 
from the same trajectories (|43)) - ([44|) : 



Co Co Co a 



i2 



-^(l+^o) 
1-Co 



2c 



1 - Co 1 - Co 



(48) 



(49) 
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Microheater location i,^ 

FIG. 3. Displacement of the collapse point relative to the 
starting point of expansion (microheater location) as deter- 
mined from numerical solution of dimensionless pump equa- 
tions (fT2|l - (fT3)l for ^ = m = 0, 71 = 72 = 1. Different lines 
correspond to different bubble strengths a — 0.1, 0.2, . . . , 1.0. 
The dashed line is the strong-asymmetry asymptote, Eq. (|47|l . 

The total post-coUapsG momentum is 
qc^L- + (1 - L) ■ i2c = 2a(l + Co) + 0{il) . (50) 

The linear ^0 term comes from the increased mass that 
the short arm picks up between return time ri^ and col- 
lision time Tc. 

Notice that the post-collapsc momentum can exceed 
twice the initial fluid momentum a, as can also be seen 
from the exact numerical solution presented in Fig. |4l 
The reason for this excess of momentum is as follows. In 
the symmetrical model, collapse dynamics is simply re- 
versal of expansion. The short arm returns to the starting 
position with a momentum equal in magnitude but oppo- 
site in sign to its initial value. If ^0 ^ Ij the return time 
is very small, so that the long arm has lost very little of 
its initial momentum. By the time the short arm returns 
to the combined momentum is almost (slightly less 
than) 2a. However, since the long arm has moved away 
from it has a little more time to accelerate before 
collision. During this extra time the short arm picks up 
more momentum than the long one loses, which results 
in an excess. Ultimately, the momentum is provided by 
the left reservoir. 



3. Complete numerical solution 

With just two dimensionless parameters, the system 
can be easily analyzed numerically to completion. The 
numerical solution is presented in Figs.[3][6l Figures[3]and 
2] show the primary and secondary effects for a number of 
bubble-strength values spanning the practically relevant 
interval 0.1 < a < 1. The most notable feature of both 




Microheater location i,^ 

FIG. 4. Post-collapse velocity as determined from numerical 
solution of dimensionless pump equations (|12p - (|13p for /3 = 
m = 0, 71 = 72 = 1. The lines correspond to different bubble 
strengths a = 0.1, 0.2, . . . , 1.0. 



sets of graphs is an optimal microheater location at which 
the effects are maximal. Note that although the optimal 
^0 for the two effects are not exactly equal, they are close, 
and share a similar trend: the optimal location is close to 
the channel's edge for weak bubbles and shifts toward the 
channel center as bubbles grow stronger. The existence 
of an optimal microheater location is a welcome feature 
for practical applications of the effect. Clearly, fabricat- 
ing microheaters right at the channel's edge would be 
challenging. The findings imply that this is unnecessary. 

Figures [5] and [6] show both inertial effects for the entire 
parameter space < ^0 < 1 and 0.01 < a < 3.0 as color 
maps. Bubble strengths a > 1 are considered impractical 
but are still of interest from the fundamental standpoint 
and are included here as such. Several features are worth 
noting, (i) For both effects, the optimal microheater lo- 
cation converges to ^0 ~ 0.35,0.65 at a > 1. (u) The 
primary effect saturates with bubble strength, although 
the a-dependence is not monotonic. (iii) The secondary 
effect does not saturate with bubble strength. It grows 
roughly linearly at large a, in agreement with Eqs. (|4ip 
and dSni). 

We have also verified that the numerical solution agrees 
with the analytical asymptotes derived earlier in this Sec- 
tion. 



C. Collapse phase. Asymmetrical model (m = 1) 

In the asymmetrical model, collapse is not simply a 
time reversal of expansion. Consider the left fiuid-vapor 
interface. After time rie, Eq. ([M)) . flow stops at ^ = 
given by Eq. After that, collapse proceeds according 
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FIG. 5. (Color online.) The complete solution of pump equa- 
tions (fT^ - (fT3t for /3 = m = 0, 7i = 72 = 1. Shown is 
the primary pumping effect: the displacement of the collision 
point relative to starting point as a function of microheater 
location and bubble strength a. The step between contour 
lines is = 0.05. 



to the equation 



FIG. 6. (Color online.) The complete solution of pump equa- 
tions HHl-ITSll for /3 = m = 0, 7i = 72 = 1. Shown is the 
secondary pumping effect: post-collapse velocity pc as a func- 
tion of microheater location and bubble strength a. The 
step between contour lines is A^' = 0.5. 



The weak-asymmetry and strong-asymmetry limit 
cases are discussed next. 



aci' + ^ef = 1, 



(51) 



with a zero initial velocity. The Bernoulli correction re- 
sults in a portion of pressure difference being spent on 
accelerating fluid within the reservoir. The implications 
are more pronounced whenever high collapse velocities 
are encountered, i.e. for strong asymmetries and strong 
bubbles. Indeed, when the mass at the turning point 
is small it accelerates quickly to a velocity of ~ 1. At 
this moment the Bernoulli correction kicks in, the pres- 
sure difference within the channel drops and acceleration 
slows down. The velocity approaches the limit value of 
\/2. New features of the asymmetrical model discussed 
below are consequences of this fact. 

In the equation of motion, Eq. (|5ip. both integrations 
are elementary with the following results 



2^^ + a ~ 



(52) 



T - Tie 



(53) 

The left edge's return velocity and time are obtained from 
here by setting = ^o: 



Tir 



Tie 



;ln 



l + \/T^ 



(54) 



(55) 



where e = exp (— q:^/2^q). Notice that return velocity 
is equal to initial velocity a/^o only in the case of weak 
bubbles, a <^1. 



Weak asymmetry 



Using ^0=2+^'' P8|) . and expanding for V ^ 1 
one obtains after some algebra the return times 



_ a Vl - ep + £o In (1 + ^1 - ep) + g^ep 



, a < 1 



2a - fa^ -( 



(56) 



(57) 



h ^^^^ + ep In (1 + ^/T^) 
2 ^/2 

+\/2epa2 |2a2 - i + 2 In (1 + Vl~^o) 



1 



(1 + ^T~^)^T^ 
lOa^ + . . . ^ a < 1 

. , a > 1 ' 



V2 ^ a 



(58) 



where ep = exp (— 2a^). The return velocities follow from 
Eq. ([54]) for smaU ip: 



Clr,2r = ± /ip - hii^ , 

ho{a) = v/2(l - ep) , 



hi{a) = 



4\/2a2eo 



(59) 
(60) 

(61) 



With return times and velocities at hand, the colli- 
sion point can be determined by linearizing trajectories 
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FIG. 7. Primary displacement in the asymmetrical pump 
model (HI])- {131) for ^ = 0, 71 = 72 = m = 1. The lines 
correspond to different bubble strengths a — 0.1, 0.2, . . . , 1.0. 



FIG. 8. Post-collapse velocity in the asymmetrical pump 
model mj-llSl) for /3 = 0, 71 = 72 = m = 1. The lines 
correspond to different bubble strengths a = 0.1, 0.2, . . . , 1.0. 



around ^o- Applying the same method as in the symmet- 
rical model the collision time and coordinate are found 
to be 



c + • V.' + Oi^p'^) 



Thus, the primary pumping effect is 



(62) 
(63) 



V2(l-6o).9V^ = | a»l. (64) 



Note that transition from the weak-bubble to the strong- 
bubble regime is not monotonic. The coefficient by if) 
passes through a minimum « —3.04 at a = 1.12 before 
converging to its large-a limit of — 1 . 

The post-collapse velocity can be derived from the 
combined mechanical momentum of the two arms at the 
time of collision. With the return times and velocities 
given above and accelerations of both arms known from 
the original dynamic equations, the velocities near the 
return times can be approximated by linear functions of 
time: 



£,[{t) = ho - hiijj + 



1 



if '2 



^0 

-1 - 



(T-c-gij), (65) 

1 t'2 
2Slr 



(r - c + ffV) • (66) 



Substituting here collision time r = c, multiplying the 
velocities by respective column lengths and (1 — ^c) 
and adding together, one obtains the total post-collapse 
momentum in the leading order in -0: 



9c 



qic + q2c = (2/10 - hlg - hi-2g)ip . (67) 



Here g, ho and hi are given by explicit expressions in 
Eqs. (|58p . (jnni) and (|6ip . Since the combined length after 



collapse is one, the post-collapse velocity is equal to the 
post-collapse momentum. The weak-bubble and strong- 
bubble limits of secondary pumping effect are 



^--^ -1^, a»l 



(68) 



Notice that unlike the symmetrical model [compare with 
Eq. (j4ip ]. here there is an optimal bubble strength {a = 
1.05), for any small microheater asymmetry ip. 



2. Strong asymmetry 

Under strong asymmetry, the microheater is close to a 
channel end, ^0 1- In this Section, condition a/^o 3> 1 
will also be assumed, meaning that bubbles cannot be 
arbitrarily weak. Physics in this regime is dominated by 
the fact that the short left arm accelerates fast during col- 
lapse and quickly reaches the limit velocity of V2. After 
that motion is uniform. Referring to general expressions 
([Ml) - ([55)) and taking into account that e is exponentially 
small, one obtains the return velocity and time 



/2 



Tlr = 



Tie 



a/2 a 



(69) 
(70) 



Thus, the return time is linear in rather than 
quadratic as in the symmetrical model [compare with 
Eq. (|42|) ]. Linearized trajectories near Tir are 



Ci(t) =^0 + V2(t. 
6(T)=eo 



Tlr) 



(71) 

(72) 



For the purposes of this Section, it is not necessary to 
include the quadratic term in the right edge's trajectory. 
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FIG. 9. (Color online.) The complete solution of pump equa- 
tions inj-HHI for /3 = 0, m = 1, 71 = 72 = 1. Shown is the 
primary pumping effect: displacement of the collision point 
relative to starting point ^0 as a function of microheater lo- 
cation ^0 and bubble strength a. The step between contour 
lines is = 0.05. 
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Microheater location 

FIG. 10. (Color online.) The complete solution of pump 
equations (P^ - (fT5)l for /3 = 0, m = 1, 71 = 72 = 1. Shown 
is the secondary pumping effect: post-collapse velocity pc as 
a function of microheater location ^0 and bubble strength a. 
The step between contour lines is A^' = 0.1. 



Equating ^1 ~ ^2 yields time and position of the collision 
point 

Co 



V2- 



a 



(73) 



(74) 



These formulas assume a < ^/2. The singularity is phys- 
ical. If initial velocity of the right edge a/{l — S,o) « a is 
close to the limit velocity of the left edge v^, it takes a 
long time for the left edge to catch up. Under this con- 
dition, the collision time and displacement are no longer 
small and higher-order terms must be taken into account. 
The last equation yields the primary pumping effect 



= - ^0 = 



(75) 



which is linear in small parameter ^o- 

The secondary pumping effect derives from the total 
mechanical momentum at collision 

qc^^cV2+{l-^c)-^ ^a + ^o{V2 + a) . (76) 
t - 4o 

Thus the limit value at ^0 ^ is a and not 2a as in 
the symmetrical model, cf. Ecj^. ([5(1)) . This is because the 
left arm has a limited velocity. The first order correction 
in ^0 is still positive which suggests the existence of an 
optimal microheater location. This conclusion is valid 
only for not very strong bubbles, a < V2. 



3. Complete numerical solution 

A full numerical solution of the asymmetric pump 
equations is presented in Figs. [TlfTUl Referring to Figs. [7] 



and[51 several differences from the symmetrical model are 
apparent. The primary effect is linear at small ^0 rather 
than quadratic. At the same time, optimal microheater 
locations and primary maxima are roughly the same as 
in the symmetrical model. In contrast, the secondary 
effect (Fig. [S]) is roughly half of the corresponding sym- 
metrical model values. Peaks are less pronounced but 
evolve with bubble strength in the same manner. Thus, 
the asymmetrical model predicts smaller net flows. 

Examination of the all-parameter solution of Figs. [5] 
and [To] reveals further differences at large bubble 
strengths. Specifically, the optimal microheater location 
shifts back toward the channel edge. In addition, both 
effects feature maxima as functions of a, implying an op- 
timal bubble strength for a given microheater location. 
Finally, the secondary effect does not grow indefinitely at 
large a, unlike in the symmetrical model, but saturates 
instead. One should note that although these qualitative 
differences between the two models are intriguing and de- 
serve additional scrutiny, they occur in the regime of very 
strong vapor bubbles which is hard to realize in practical 
devices. 



Summarizing Section IIIIl one concludes that both 
models studied describe robust inertial pumping. Rel- 
ative simplicity of the models has allowed many results 
to be derived analytically, which has elucidated physics 
behind the effect. There are enough differences in predic- 
tions that a comparison with CFD calculations is likely 
to suggest the most realistic model and the underlying 
channel-reservoir boundary condition. Insights gained 
from this analysis will be helpful in designing practical 
inertial pumps. 
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FIG. 11. The complete solution of pump equations (|12p - (|13p 
for the symmetrical model m — and 71 =72 ~ 1- The 
bubble strength is a = 0.5. Shown is the primary pumping 
effect as a function of microheater location ^0 for different 
friction parameters /3 — 0.0,0.5, 1.0, 10.0. Trace with the 
largest amplitude corresponds to 13 — 0.0. 



FIG. 12. Same as Fig. [TT] but the secondary pumping effect 
is shown. 



ertial pumping. 



IV. 



EFFECTS OF VISCOSITY 



So far, viscous forces have been neglected. It is intu- 
itively clear that for a mechanism that relies on fluid iner- 
tia viscosity should have detrimental effect. Within the 
one-dimensional dynamical model viscosity is included 
via dimensionless parameter /3 defined in Eq. p^ . j3 is 
directly proportional to bulk viscosity of the participat- 
ing fluid, cf. Eq. To get a sense of the physically 
relevant /? interval, a fluid with 77 = 1.3 mPa-s in a mi- 
crochannel 200 //m long and 20 fiiii wide and tall would 
have /3 w 2. 

A non-zero 13 introduces additional non-linearities into 
the equation of motion. Although some useful results 
can still be obtained analytically (particularly, in the 
post-collapse phase, see Section |V|, final pumping effects 
are difficult to derive. We therefore resort to numerical 
analysis. Results are presented in Figs. [11] and [T^] for 
the symmetrical model and in Figs. [13] and [14] for the 
asymmetrical model. In selecting a representative bub- 
ble strength, we chose a = 0.5 which roughly corresponds 
to vapor bubble pressure p^, = 8 — 10 atm, (cf. caption 
to Table [J). 

As expected, finite viscosity systematically reduces 
both inertial effects. However, reduction is not uniform 
across all microheater locations. Reductions are stronger 
for 0.1 < ^0 < 0.9 but are much less for the locations 
close to the channel ends, as can be observed in the 
plots. For the primary effect. Figs. [TT] and [13] such a 
non-uniform change results in a systematic shift of the 
optimal location closer to the end of the channel. For the 
secondary effect, the maximum is eliminated altogether. 
This behavior is another indication of complexity of in- 



V. POST-COLLAPSE PHASE 

For realistic conditions, most of the net fiow happens 
in the post-collapse phase when fluid moves by inertia 
against friction forces. It is therefore important to know 
details of the post-collapse flow. The equation of motion 
was derived in Section HlDl Eq. (|^ . Since this equation 
is simpler than its counterparts of the expansion-collapse 
cycle, it is possible to obtain analytical solutions even for 
non-zero frictions /3 and pressure heads 71 — 72. In this 
Section, the most common case of zero pressure head is 
solved. The more complex 71 ^ 72 case will be considered 
in a separate study. 

For the symmetrical model, m = 0, the equation of 
motion reduces to 



with the obvious solution 



e'(r)=77,e-^(^-^=) 
e(r)^ec + |[l- 



,-/3(T-r,) 



(77) 

(78) 
(79) 



Flow slows down according to a simple exponential law. 
The total displacement of the post-collapse phase is rjc/ (3. 
It diverges at zero viscosity. 

In the asymmetrical model, m = 1 , the Bernoulli term 
changes sign depending on flow direction. However, in 
the case of zero pressure head, flow direction does not 
change with time. It is sufficient therefore to consider 
only one, for example, positive flow direction, ^'(r) > 0. 
The dynamic equation then reads 



(80) 
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FIG. 13. The complete solution of pump equations (|12p - (|13p 
for the asymmetrical model m = 1 and 71 = 72 = 1- The 
bubble strength is a = 0.5. Shown is the primary pumping 
effect as a function of microheater location ^0 for different 
friction parameters /3 — 0.0,0.5, 1.0, 10.0. Trace with the 
largest amplitude corresponds to /3 — 0.0. 



Both integrations are elementary yielding 



{?],, + 2/3) e'3(— c) - rj. 



e(r) =C, + 21n 



r]c_ 
2/3 



1 



(81) 

(82) 



The total displacement of the post-collapse flow is finite 



(83) 



A^(oo) = 21n(^ + l 



In the zero- viscosity limit, velocity decays to zero as 
2r]c/[2 + 'r]o{T — Tc)], but the overall displacement diverges 
logarithmically. 



VI. 



SUMMARY 



A micropump is required for any active fluid-handling 
system. Inertial micropumps do not contain moving 
parts and can be made in large quantities by batch fab- 
rication processes. As such, they are an excellent can- 
didate for the universal integrated pump of chip-scale 
fluidics. Physics behind pump operation is defined by 
a subtle balance between pressure, viscous, and inertial 
forces. A high-pressure vapor bubble generated by a mi- 
croheater expels fluid from the channel to reservoirs, but 
after that flow reverses under a negative pressure differ- 
ence. Because of a smaller mass and inertia the shorter 
arm reverses flow direction earlier and then has more time 
to accelerate during inflow. By the time of collision, the 
shorter arm acquires larger velocity and momentum than 
the longer arm. The two columns collide at a point that 



FIG. 14. Same as Fig. [13] but the secondary pumping effect 
is shown. 



is shifted from the initial point of expansion, which con- 
stitutes the primary pumping effect. A non-zero total 
momentum ensures post-collapse flow from the short to- 
ward the long side of the channel, which is the secondary 
pumping effect. Total net flow is the sum of the primary 
and secondary contributions. 

In this paper, inertial pumping action has been stud- 
ied within a simplifled onc-dimcnsional dynamic model. 
Transverse motion within the channel has been neglected 
and entire dynamics reduced to that of a fluid-vapor in- 
terface treated as a mathematical point. Despite these 
simplifications, the one-dimensional model captures main 
features correctly, while using scaled units reduces the 
number of independent parameters. The model becomes 
an efficient tool of analyzing different pumping regimes. 
A major challenge lies in understanding pressure at the 
channel-reservoir interface. Transient nature of the flow, 
vortex formation in the reservoir, non-uniform velocity 
across the channel cross-section during flow reversal, and 
other factors make selection of a single boundary condi- 
tion non-trivial. Until full understanding is achieved via 
additional experimental and numerical work, the current 
approach is to assume a physically reasonable boundary 
condition and analyze its consequences. In this paper, 
two possible choices have been studied in detail. In the 
simplest "symmetrical" model the interface pressure is 
assumed constant and set to be equal to the bulk reservoir 
pressure. In the more complex "asymmetrical" model, 
during inflow the interface pressure is reduced from the 
bulk value by a Bernoulli-type correction. 

A key finding of the paper is that both models con- 
tain pumping effects and their properties are similar for 
weak-to- intermediate bubble strengths (which is the most 
realistic situation). In particular, the two models predict 
an optimal microheater location in the 0.2 — 0.3 range 
for both the primary and secondary effects. A non-zero 
viscosity dampens net flow in both models in a similar 
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fashion. The predictions begin to diverge for very strong 
bubbles a > 1, as can be seen by comparing Figs. [SJ [H] 
with Figs, m [ini The approach adopted in this paper 
can be extended in a number of ways including sequen- 
tial firing, non-zero pressure heads and branched chan- 
nels. These systems will be investigated in subsequent 
publications. 
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Appendix A: Derivation of Eq. (f6)) . 

For Poiseuillc flow in a cylindrical channel of radius R, 
the velocity profile and the total fiow are given by 



W 



8^^ ' 



(Al) 



(A2) 



where 77 is the bulk viscosity of the fluid. Connection to 
the one-dimensional model studied in the main text is 
done by identifying cross-section average velocity W/A 
with fluid velocity x of the one-dimensional model. Ex- 
pressing from here the pressure difference 



Ap 



8x 



(A3) 



and substituting into Eq. (|Aip . the velocity profile is 
rewritten as 



v^{r) = 2x 



(A4) 



The only non-zero component of the stress tensor is: 

The total viscous force acting on a cylinder surface of 
length X is 



\r=R 



■ 2'nRx = ~8TTr]xx . 



(A6) 



Comparison with the definition of the viscous force term, 
Eqs. (HI)-®, yields the relationship ©. 

Notice that the dimensionless conversion factor Stt « 
25 is unusually large. For example a 30% solution of 
glycerol in water at 45 °C with 77 = 1.3 mPa-s, would 
have an effective friction coefficient of k w 32 mPa-s. In 
[l| we reported a value of k = 28 mPa-s extracted from 
fits to full-scale CFD simulations, which is close to the 
above estimate. 
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